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Abstract  A  n^w  algorithm  is  presented  for  computing  vertices  of  a  simplicia!  ^|-/ 
(nangulation  of  the  p-dimensional  solution  manifold  of  a  parametrized  equation 
F(x)  =  0,  where  F  is  a  nonlinear  mapping  from  Rn  to  R™,  p=n-m>1 .  An  essential  part  of 
the  method  is  a  constructive  algorithm  for  computing  moving  frames  on  the  manifold; 
that  is,  of  ortho  normal  bases  of  the  tangent  spaces  that  vary  smoothly  with  their  points 
of  contact  The  triangulation  algorithm  uses  these  bases,  together  with  a  chord  form  of 
the  Gauss- Nev»ton  process  as  corrector,  to  compute  the  desired  vertices.  The  Jacobian 
matrix  of  the  mapping  is  not  required  at  all  the  vertices  but  only  at  the  centers  of  certain 
localTriangulation  patches".  Several  numerical  examples  show  that  the  method  is  very 
efficient  in  computing  triangulations,  even  around  singularities  such  as  limt  points  and 
bifurcation  points  This  opens  up  new  possbilties  for  determining  the  form  and  special 
features  of  such  solution  mat  .ifolds. 
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1  Introduction 


Parameter-dependent  nonlinear  equations 


F(2,A)  =  0,  (1.1) 

involving  a  state  variable  z  and  a  parameter  vector  A,  arise  in  many  applications 
Under  natural  cond lions  on  F  and  the  relevant  spaces  the  set  of  solutions  (z.  A)  of 
(1.1)  constitutes  a  differentiable  manifold  in  the  product  of  the  state  and  parameter 
space,  and  the  dimension  of  this  manfold  equals  the  parameter  dimension. 

In  most  practical  applications  interest  centers  not  so  much  on  computing  a  few 
solutions  of  (1.1),  but  rather  on  determining  the  form  and  special  features  of  the 
solution  manfold.  For  instance,  f  (1 .1)  represents  an  equilibrium  problem,  then  we 
may  wish  to  determine  the  bfurcadon  diagram  or  the  boundaries  of  the  stabilly  regions 
on  the  manfold  But,  as  1  turns  out  all  standard  computational  methods  for  such  an 
analysis  require  us  to  construct  a  picture  of  a  p-dimensional  manfold  from  information 
along  one-dimensional  paths.  In  fact,  all  these  methods  belong  to  the  family  of 
continuation  processes  for  which  the  dimension  of  the  parameter  space  always  has  to 
equal  one  Thus,  before  such  a  process  can  be  applied,  any  problem  with  a  larger 
parameter  dimension  must  be  reduced  to  some  form  involving  only  a  scalar  valued 
parameter  and,  geometrically,  such  a  reduction  is  equivalent  wth  a  restriction  to  some 
path  on  the  solution  manfold  of  the  original  equation  A  continuation  method  then 
computes  a  sequence  of  points  along  such  a  path.  For  example,  in  structural 
engineering  the  parameter  A  often  characterizes  a  vector  of  load  components  in  which 
case  1  has  become  customary  to  fix  a  linear  combination  of  these  components 


speciying  a  particular  load  direction  The  resulting  reduced  equation  then  involves 
only  the  load  intensity  as  a  one-dimensional  parameter  variable  and  the  standard 
incremental'  methods  generate  points  along  this  load  path. 

In  general,  it  is  not  easy  to  develop  a  good  picture  of  a  p-dimensional  ma  nifold  from 
information  along  one-dimensional  paths  Thus  it  is  not  surprising  that  there  is  growing 
interest  in  computational  methods  which  generate  multi-dimensional  -grids  of  solution 
points  covering  an  entire  segment  of  the  manifold.  Up  to  now  the  only  method  for 
computing  such  multi-dimensional  grids  appears  to  be  that  of  E.L.AIIgower  and 
P  H  Schmidt  (2).  It  utilizes  a  simplictal  continuation  algorithm  to  "triangulate",  by  means 
of  p-smnplices,  a  portion  of  a  p-dimensional  manifold  defined  by  an  equation  of  the 
form  (11) 

Here  we  present  a  different  method  for  computing  vertices  of  such  a  triangulation  of 
segments  of  the  p-dimensional  solution  manifold  of  an  equation  (1  1)  An  essential 
part  is  a  constructive  algorithm  for  competing  ortho  norma.  I  moving  frames  on  the 
manifold  in  the  sense  first  considered  by  E.Cartan,  that  is,  of  ortho  norma  I  bases  of  the 
tangent  manifolds  that  vary  continuously  wth  their  points  of  contact  (see  e  g  ,  (1 7] ) 

The  resulting  triangulation  algorithm  uses  these  bases  and  a  predictor-corrector 
approach  to  compute  the  desred  grid  points  Its  has  ma ny  similarities  with  the 
continuation  methods  including  a  comparable  computational  complexity  In  particular, 
trie  Jacob®  n  matrix  of  the  mapping  is  not  required  at  all  the  points  For  example,  on  a 
two-dime  ns  io  rial  manifold  the  computation  of  a  typical  triangulation  pattern  with  I  H 
triangles  requires  only  1 9  Jacobian  decompositions 

Alter  summarizing  some  basic  concepts,  in  Section  2  we  introduce  the  moving  frame 
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algorithm  in  Section  3  Then  Section  4  outlines  the  general  Iriangulation  method  arid 
Section  5  presents  several  numerical  examples.  Finally,  we  end  with  an  outlook  on  the 
utilization  of  the  computed  tnangulations  for  the  determination  of  specific  features  of 
the  manifold 
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2  Basic  Concepts 
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Throughout  this  article,  let 

F:S  -e  R™,  S  open  in  Rn,  p  =  n-m>1,  (2.1) 

be  a  differentiable  mapping  of  class  Cr,r>  1,  on  the  subsets.  As  usual,apointx  s  S 
is  called  regular  i  the  first  derivative,  DF(x),  of  F  has  full  rank  m  and  hence  maps  Rn 
onto  Rm  We  consider  the  equation 

F(x)  =  0,xeS,  (2.2) 

and  assume  that  is  regular  solution  set 

M  =  { x  e  S ,  F(x) «  0,  x  regular)  (2.3) 

is  non-empty,  t  is  well-known  that  M  is  a  p-dimensional  CT -manifold  in  Rn  wthout 
boundary  (see,  e  g.,  [17]  or  [15]). 

The  tangent  space  T^  at  any  point  x  e  M  may  be  identified  with  the  kernel  of  the 
Jacobian  DF(x),  that  is, 

TyM  =  ker  DF(x)  =  {  u  e  Rn ,  DF(x)u  =  0).  (2.4) 


The  norma!  space  NX,M  at  x  e  M  is  the  orthogonal  complement  of  the  tangent  space 


under  the  natiral  Euclidean  inner  product  on  Rn,  that  is, 

NxM  =  (TxM)i  =  (kerDF(x))1  =  rgeDF(x)T  (2.5) 

Since  DF(x)  has  maximal  rank  in  some  open  subset  Sgof  S  containing  M,  the 
mapping 

xeS0  DF(x)T[DF(x)DF(x)T]'1DF(x)  e  L(Rn)  (2  6) 

from  S0  into  the  space  L^R^ofall  linear  isomorphisms  on  Rn  is  of  class  Cr*1  on  Sq. 
Hence,  the  orthogonal  projection 

P  M  L(Rn);  P(x)  =  ln  -  DF(x)t  (DF(x)DF(x)Tr1  DF(x),  x  e  M  (2.7) 

of  Rn  onto  T^M  is  a  mapping  of  class  C*"'1  on  the  manifold  M,  (here  ln  denotes  the 
identity  on  Rn). 

For  the  computation  we  require  local  coordinate  systems  on  M.  Any  p-dimensiona  I 
subspace  Tof  Rn  induces  a  local  coordinate  system  of  M  at  any  point  x  e  M  where 

T  n  NXM  =  {0},  (2.8) 


In  fact,  if  (2  8)  holds  forx  e  M  then  there  exist  open  neighborhoods  Y-|  andY20fthe 
origins  of  T  and  Rn,  respectively,  as  well  as  a  unique  Cr'1  function  w  V1  ->  T1  with 
w(0)  =  0,  such  that 


MnY2  =  {ye  Rn;  y  =  x  +  t*v<t),  teY^,  (2.9) 

(see,  e  g.,  [9J  or  [1 5J).  In  other  words,  in  the  local  coordinate  system  induced  by  T ,  the 
point  y  *  x+t+w(t)  of  M  has  the  coordinate  t  e  T. 

A  point  x  e  M  where  (2.8)  holds  is  a  non-singular  point  with  respect  to  the  given 
coordinate  space  T,  else  we  call  x  a  singular  point  Clearly,  at  any  point  x  e  Mthe 
tangent  space  can  be  chosen  as  coordinate  space  and  x  is  non-singular  with 

respect  to  it.  In  most  applications,  a  "natural"  parameter  space  A  is  given,  as  indicated 
by  the  form  of  the  equation  (1 .1),  and  the  orthogonal  subspace  Z  =  a-*-  is  the  state 
space  Then  interest  centers  on  determining  the  singular  points  wfth  respect  to  the 
space  A  These  are  the  so-called  foldpoints  on  M  where  the  tangent  space  has  a 
non-zero  intersection  with  the  state  space  Z  and  the  parameter  space  A  can  no  longer 
be  used  as  a  local  coordinate  space.  These  are  also  the  points  where,  for  example,  in 
equilibrium  problems  a  change  in  the  stabiliy  behavior  of  the  physical  system  under 
study  may  be  expected. 


Numerically,  the  mapping  w of  (2. 9)  can  be  rnplemented  in  various  ways  A  smple 
approach  is  based  on  a  chord  form  of  the  Gauss- Neviton  method.  At  the  given  point 
x  e  M  we  compute  the  QR-factorization 


DF(x)T 


(2.10) 


of  the  transposed  Jacobian  DF(x)T  involving  the  n  x  n  orthogonal  matrix  Q  and  m  x  m 
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nonsingular,  ipper  triangular  matrix  R.  Then,  starting  from  any  point  y  in  a  suitable 
neighborhood  of  x  in  x+T^,  we  may  apply  the  process 


1)  Sety°  =  y; 

2)  for  k  =  0,1,...  until  convergence 
2a)  solve  RTz  =  F^)  for  z  e  RP; 

2b)  compute  the  next  iterate  y*+1  =  yk-Q  (z,0)T; 

With  (2.1 0)  this  is  readily  rewritten  in  the  form 
DKxXy^’-yk)  +  F(/)=  0  ,  y^-y*  e  [kerDF^j1,  k=0,1, ... 


(211) 


(2.12) 


which  shows  that  yk+1-yk  e  tyvi  Thus  we  have  y^e  y+N^  for  all  k>0, 
whencethe  limt  pointy*  -  if  l  exists-  is  the  unique  point  in  the  intersection  of  M  and 
y+hyvl  which,  in  the  notation  of  (2.9),  can  be  written  as  y*=x+t+w(t),  t=y-x. 

The  convergence  theory  of  Gauss-Newton  processes  is  well  understood.  Earlier 
studies  of  these  methods  considered  applications  to  least  squares  problems  and 
hence  assumed  that  F  maps  Rn  into  R™  where  n  <  m.  A  local  convergence  result  which 
covers  our  case  n  >  m  may  be  found  in  (8,  Theorem  4],  Another  simple  proof  for  the 
method  in  the  form  (21 2)  also  follows  along  the  lines  of  the  convergence  proof  for 
singular  chord  methods  given  in  [14]  These  results  guara  ntee  the  validity  of  the 
following  theorem 

Theorem  1  Under  the  stated  assumptions  about  the  mapping  F  there  exists  for  any 
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3  A  Moving  Frame  Algorithm 

As  usual,  a  vector  field  of  class  Cs,  s  s  r,  on  an  open  subset  Mq  of  our  manifold  M  is 
a  Cs  function  u:  Mg  -» TM  from  Mg  into  the  tangent  bundle  TM  such  that  u(x)  e  T^for 
all  x  e  M0  Amoving  frame  of  class  Cs  on  Mg  is  a  mapping  which  associates  with 
each  x  g  Mg  a  frame  (i  e.,ordered  basis)  {u1, uP}  of  such  the  functions 
u'  Mg  ->TM,  i*1,.  p,  form  p  vector  fields  of  class  CsonMg.  When  such  a  moving 
frame  exists  on  Mg  then  the  sub-manifold  Mg  is  said  to  be  parallelizable  We  will 
consider  only  orthonormal  moving  frames;  that  is,  frames  for  which  the  basis  vectors 
are  orthonormal.  (Fora  discussion  of  these  concepts  see,  e  g.,  [17]  ) 

Clearly,  the  problem  of  computing  an  orthonormal  basis  of  the  tangent  space  TXM  of 
Mat  a  given  point  x  g  M  is  equivalent  wih  the  construction  of  an  n  x  pmatrix  U  with 
orthonormal  columns  for  which 

DF(x)U  =  0  (31) 

There  are  many  techniques  for  computing  such  a  matrix.  A  well-known  procedure  is 
provided  by  the  QR-decomposition  (2.1 0).  In  (act,  if  the  matrix  Q  is  partioned  in  the  form 
Q  =  (Q-|,  Q2)  where  Q-j  has  m  columns  then  we  may  use  U  =  Q2  as  the  desired  basis. 
Various  other  techniques  for  computing  U  have  been  proposed.  We  mention,  in 
particular,  the  methods  in  [5]  and  [6]  aimed  at  producing  sparse  bases  U 

An  algorithm  for  constructing  a  moving  frame  of  class  Cs  on  some  open  subset  Mg 
of  M  has  to  generate  a  basis  matrix  U  =  U(x)  for  each  x  e  Mg  in  such  a  way  that  the 
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mapping  U:  Mg  -»  L(RP,  R”)  is  of  class  Cs  As  T.F.Coleman  and  DC.  Sorensen  [7] 
have  noted,  the  approach  based  on  the  QR-decompostoon  does  not  give 
continuously  varying  matrices  U(x).  This  observation  extends  to  other  a  Igor  thms  of  a 
similar  nature;  in  fact,  I  relates  directly  to  the  corresponding  problems  of  computing 
eigenvectors  associated  wth  a  mutple  eigenvalue.  Three  remedies  are  proposed  in 
[7j,  but  they  concern  only  the  construction  of  a  liml  Ug  of  a  sequence  of  bases  U(x) 
when  x  tends  to  x0. 

For  our  construction  of  a  moving  frame  we  restrict  attention  to  open  subsets  Mg  of  M 
where  a  given  p-dime  ns  tonal  subspace  T  of  Rn  induces  a  local  coordinate  system,  that 
is,  where  (2.8)  holds  for  all  x  of  Mg.  Wth  a  mild  abuse  of  language,  let  Tgbeapx  p 
matrix  with  orthonormal  columns  which  span  the  coordinate  space  T.  Moreover, 
assume  that  we  have  picked  some  method  for  computing  at  any  point  xofM0annx  p 
matrix  U(x)  wth  orthonormal  columns  that  span  the  tangent  space  TXM  at  x  Of  course, 
U  is  not  expected  to  depend  continuously  on  x.  For  instance,  we  may  use  the  matrices 
produced  by  the  QR-decompostion  technique.  For  any  orthogonal  pxp  matrix  Q  = 
Q(x)  the  matrix  U(x)Q  is  another  ortho  norma  I  basis  of  T^  and  our  am  is  to  construct 
matrices  Q(x)  such  that  the  ’rotated"  bases  U(x)Q(x)  depend  continuously  on  x  for  all  x 
in  M0 

The  normalization  (T0)TT0  =  Ip  suggests  that  we  choose  the  orthogonal  matrix  0 
so  that  (U(x)G)tT0  approximates  the  identity  Ip.  Various  norms  may  be  used  for  this, 
an  advantageous  choice  is  the  Frobemus  norm  ||  A  |jp  =  [  tr  (A^A)  ] 1/2  The  resulting 
optimization  problem 


||(U(x)Q)TT0-ip)lF  =  mm  ,  subject  to  QTQ  =  lp 


(3  2) 


s  a  case  of  the  orthogonal  Procustes  problem.  As  discussed  In  (1 2],  the  following 
algorithm  solves  (3.2): 

(I)  U0  :=  UWTT<>, 

(  2)  compute  the  singular  value  decomposfoon  (3  3) 

atU0B  =  Z  arid  save  A  and  B. 

(3)Q  =ABt 


For  our  purposes  the  essential  fact  is  now  the  content  of  the  following  theorem 

Theorem  2  Let  be  an  open  subset  of  Mon  which  the  given  p-dimensional 
sifcspace  Tof  Rn  induces  a  local  coordinate  system  For  any  x  e  M^  let  U(x)  be  an 
ortho  normal  basis  matrix  of  Tj<M  and  compute  the  orthogonal  matrix  Q  =  Q(x)  of  (3  3) 
Then  the  mapping  x  e  -»  U(x)Q(x)  e  L^RPj}11)  is  of  class  CT'1  on  M0  and  defines 
an  orthonormal  moving  frame  on  Mq 

Proof  Evidently  ,  UqZ  =  0  implies  that  the  tangent  vector  U(x)z  e  T^M  must  be 
orthogonal  to  the  subspace  T  of  Rn  spanned  by  the  columns  of  T0,  and,  and  that 
U(*)z  eW  By  construction  of  M0  this  cannot  happen  for  x  e  M0  unless  z  =  0  In  other 
words,  for  *  e  M0  the  matrices  U0  and  I  arising  in  (3  3)  are  non-singular  Now 


U0  =  AZBT  =  ABT(6TiB)  =  OH,  H  =  BTIB. 


ts  the  polar  decomposition  of  U0  and  it  follows  that 
H  -  l(U0)TU0l 1,2  -  [(T0)TU(*)U(x)TT0|'12 
ts  non-singular,  whence 

0  -  (U0)TT0  |(T0)T^)U(x)TT0r''2 

Evidently,  U(x)U(x)T=  P(x)  is  the  orthogonal  projection  (2.7)  from  Rn  onto  T„M  Thus 
we  see  that 

U(X)Q  =  P(x)T0  t(T0)Tp(x)T0r  "2  (3.4) 

and,  since  P  was  afready  shown  to  be  of  class  Cf' 1  on  M,  the  resut  follows. 

Oir  overall  moving  frame  algorithm  on  Mq  now  consists  of  the  following  three  steps 

(1)  Given  x  e  compute  the  basis  matrix  U(x)  of  T^; 

(2)  compute  the  orthogonal  matrix  Q  by  (3.3),  (3.5) 

(3)  form  the  desired  basis  matrix  U(x)Q. 

f  the  QR-factonzation  is  used  in  step  ( 1 ),  then  the  order  of  the  number  of  floating 
point  operations  regured  is  as  follows: 

Computation  of  U(x)  0(nm2) 

MuRplication  U(x)TT0  0(np2) 


Singular  value  decomposition  0(pJ) 

Formation  of  the  product  U(x)Q  0(n+p)p2) 

Thus,  when  the  dimension  p  of  the  manifold  is  small  in  comparison  wth  the  space- 
dimension  n,  as  is  typical  in  applications,  then  the  principal  cost  is  related  to  the 
QR-factorization  of  DF(x)T  and  involves  about  (2  /3)n3  operations.  This  is  indeed 
analogous  to  the  complexly  of  a  standard  continuation  process. 

For  the  practical  rnplementation  I  is  certainty  desrable  to  choose  the  basis  matrix 
T0  of  our  reference  coordinate  space  T  in  Rn  as  simply  as  possble.  In  particular,  I  is 
very  advantageous  to  define  T  as  a  subspace  spanned  by  p  appropriate  natural  basis 
vectors  e1  ,...,en  of  Rn  Then  T0  can  be  taken  as  a  matrix  wth  columns  e'  with  certain 
distinct  indices  hiy  Ui^n,  For  the  choice  of  these  indices,  recall  that  for  any  x 

of  M  and  given  vector  a  e  Rn  the  princpal  angle  <x  s  [0,ti/2]  between  TjjM  and 
span{a}  is  defined  by 

cos(<x)  *  max{  uTa ;  u  e  TJA,  |(u||2  =  1 } 

Evidently,  if  a  is  one  of  the  global  basis  vectors  of  Rn,  then 

cos(ai)=||U(x)Tel||2,  M,...n;  (3.6) 

that  is,  the  Euclidean  norm  of  the  i-th  row  of  U(x)  is  the  cosine  of  the  principal  angle 
between  the  tangent  space  TXM  and  the  i-th  coordinate  line  span{e'}  of  Rn.  Since  the 
Euclidean  norm  is  invariant  under  orthogonal  transformations,  I  is  obvious  that  the 
principal  angle  does  not  depend  on  the  particular  basis  U(x)  of  TXM. 


This  suggests  the  desired  selection  of  T0.  We  intialize  otr  moving  frame  algorithm 
at  some  reference  point  x*  of  M  and  compute  a  basis  matrix  U(x*)  of  the  tangent  space 
of  M  at  this  point  This  allots  for  a  straightforivard  calculation  of  the  principal  direction 
cosines  tj  *  cos(ctj)  of  (3.6)  which  we  order  in  descending  order  of  size.  Let  i^.Jp  be 
the  indices  of  the  p  largest  of  these  Tj  (wth  ties  broken, say,  by  lexicographic  ordering), 
then  the  corresponding  natural  basis  vectors  of  Rn  span  our  reference  coordinate 
space  T  and  form  the  columns  of  the  basis  matrix  T0.  S  ince  Ufx*)  has  rank  p,  none  of 
the  selected  coordinate  directions  can  be  orthogonal  to  the  tangent  space  Hence,  as 
requred,x*  is  a  non-singular  point  wth  respect  to  T,  and  the  subset  fv^of  Theorem  1 
contains  an  open  neighborhood  of  x*onM.  Geometrically  the  constructed  subspace  T 
is  close  to  the  tangent  space  of  M  at  x*  in  the  sense  of  the  above  maximization  of  the 
direction  cosines  (3.6).  In  fact,  our  construction  is  analogous  to  the  local  coordinate 
selection  used  in  the  continuation  program  PITCON  (see  [1 5]  or  [1 6]).  Obviously,  this 
choice  of  T0  also  has  the  advantage  that  the  computation  of  the  matrix  U0  in  step  (1 )  of 
the  algorithm  (3.3)  simply  becomes  an  extraction  of  p  of  the  rows  of  U(x). 

A  frequently  occurring  case  are  manifolds  with  dimension  p  =  2.  For  them  we  can 
compute,  up  to  signs,  a  direct  representation  of  the  orthogonal  matrix  G  of  the  algorithm 
(3.3).  In  fact,  sippose  that 

fa  bl 


is  the  matrix  in  step  (1)  of  (3  3)  Then  a  straightforward  calculation  shov»s  that 


t=  ±(a  +  dys ,  o  =  ± (b - cys ,  6  “[(a+d)*  +  (b-c)*]''f 


The  signs  are  readily  chosen  by  comparing  the  directions  of  the  computed  frame 
with  the  directions  of  the  natural  basis  vectors  used  in  T0.  This  leads  also  to  a  dir  ect 
proof  of  the  theorem  in  this  special  case. 


»  -I 

t  I 

4  A  Inangujation  Algorithm 

The  resuts  of  the  prey  to  us  section  are  now  used  to  generate  the  desred 
trangulation  on  a  subset  of  our  manfold  M.  The  basic  idea  is  as  follows:  We  introduce 
a  reference  triangulation  on  RP  and  use  the  bases  produced  by  the  moving  frame 
algorithm  to  map  segments  of  it  onto  the  spaces  x+TjjM  corresponding  to  appropriate 
points  xon  M  Then  the  Gauss- Newton  process  (2  11)  is  applied  to  “project  these 
trangulations  from  x+TxM  ontoM 

The  reference  triangulation  ,  of  course,  is  any  covering  of  RP  by  a  locally  finite 
collection  of  p-smplices  such  that  any  two  of  these  smplices  intersect  ether  in  a 
common  face  or  not  at  all  The  Iterature  on  this  topic  is  large  and  we  refer  here  only  to 
the  discussion  of  various  numerically  efficient  triangulations  in  (18)  Our  algorithm  does 
not  place  any  particular  restrictions  on  the  choice  of  this  triangulation  except  that  we 
should  be  guided  by  considerations  of  computational  simplicty.  Let  I  be  the  collection 
of  smplices  of  this  triangulation 

Mosttriangulations  used  in  smplicial  continuation  studies  are  generated  by  pivoting 
rules  A  simple  such  rule  is  pivoting  by  reflection  For  any  index  je  {1,2,..  ,p)  set  j+  = 
j*1  and  j_  =  j-1  with  the  provision  that  j*=  1  fj-pandj.«pfj«1.  Then,  fora 
given  p-smplex  o  =  [y°,y  ’ .  ,yP]  in  RP,  pivoting  by  reflection  of  the  v^itex  yJ  is  defined 
as  the  replacement  of  a  by  thesmplex  (y°,  ,yf'1,yJ*«-yl'-yl.yM.  yp)  fo0  tsa 
given  reference  smpiex  m  RP.  then  by  repeated  application  of  this  procedure  a 
triangulation  of  RP  can  be  generated  (see,  eg  }1J) 


A  frequently  used  example  is  the  so-called  Kuhn  tria  regulation  which  is  generated  by 
repeated  pivoting  by  reflection  starting  with  the  smplex 

o0  =(  0,e*,e*  *e2, . .^e1  ♦  e2  ♦  ...•*■  eP  ]. 

In  the  olten  occuring  case  p  =  2,  we  can  also  use  triangulations  of  R2  by  equilateral 
(riangles  generated  by  pivoting  by  reflection  beginning  wth  the  2-simplex 

oo-[0,e1,0.5(e1  W3e2)]  (4.1) 

Let  ^  denote  a  given  vertex  of  the  triangulation  inFtf>and5>0a  fixed  steplength. 
Then  for  any  point  x  e  M  where  a  basis  matrix  U  of  TXM  is  known,  the  mapping 

A  :  RP  x+TjfM ,  Arj  =  x  ♦  6U(rj-^),  q  e  RP  (4.2) 

transfers  i  onto  x+T^  As  in  Theorem  1 ,  let  Y(x)  denote  the  local  convergence 
domain  in  x+TjjM  of  the  Gauss-Neviton  process  (2.11)  If  q  is  a  vertex  of  J  for  which 
An  belongs  to  Y(x),  then  (2.11)  can  be  used  to  map  An  into  a  point  y  e  M.  The  set 
r(^x,U,S)  of  vertices  of  2  that  can  be  mapped  onto  M  in  this  way  shall  be  called  the 
'patch'  corresponding  to  the  information  ^,x.U,6 

An  idealized  version  of  our  algorithm  can  now  be  formulated  as  follows. 

( 1 )  Select  a  reference  vertex  ^  of  2 ; 

(2)  select  a  reference  point  x  =  x*  of  M  ; 


(3)  initialize  the  moving  frame  algorithm  at  x*  and  let  M0  be  the 
subset  where,  by  Theorem  2,  this  algorithm  applies; 

(4)  mark  the  vertex  ^  as  used;  (4 . 3) 

(5)  while  x  belongs  to 

(5a)  mark  $  as  a  ’center , 

(5b)  compute  the  frame  U  =  U(x)  by  the  moving  frame  algorithm, 

(5c)  select  all  vertices  of  the  patch  r(^U, 5)  which  have 
not  yet  been  marked  ’used" ; 

(5d)  use  (4.2)  to  map  these  vertices  onto  x+TyM  and  mark  them  ’used"; 

(5d)  use  the  Gauss-Neviton  process  to  project  the  result ng  points 
from  x+T^  onto  M ; 

(5e)  choose  a  "used"  vertex  £  of  I  not  marked  a  ‘center 
and  let  x  be  Is  computed  image  on  M; 

The  points  computed  on  M  inherithe  connectiYiy  pattern  of  the  original  simplices  of 
I  which  in  turn  induces  a  simplicial  approximation  M2  of  a  segment  of  M  in  Rn  The 
algonthm  is  still  "ideal”  in  nature  since,  in  practice,  the  sets  Mq  and  r(^x,U,S)  are  not 
known  explictly  Without  this  knowledge  the  computation  may  halt  when  the  iteration 
in  step  (5d)  fails  to  converge,  that  is,  when  we  encounter  a  point  in  the  affine  space 
«+TxMofoneofthe  centers  x  which  does  not  be  long  to  the  neighborhood  V(x) 
spectfied  in  Theorem  2  .  A  second  poss bilfty  for  failure  arises  in  the  execution  of  the 
moving  frame  algorithm  in  step  (5b)  when  the  selected  center  x  does  not  belong  to  M0 


In  order  to  make  the  algonthm  practical,  we  replace  the  "ideal"  patches  r(^,xU,6)  in 
step  (5c)  by  "standardized"  patches  r0(^)  The  definition  of  these  patches  depends  on 


the  specific  reference  triangulabon  in  f#\  As  an  example,  consider  the  earlier 
mentioned  tnangulation  on  R2  consisting  of  equilateral  triangles  produced  by  pivoting 
by  reflection  from  the  triangle  (41)  Then  the  standarized  patch  for  the  center  point 
(0,0)  in  Figure  1  is  the  hatched,  star-shaped  region  and  for  any  other  vertex  it  is 
obtained  by  obvious  translation  With  this  standard  patch  the  progress  of  the  algorithm 
is  easily  follloYved  in  Figure  1 .  There,  at  each  vertex,  the  second  of  the  two  integers  is  a 


counter  and  the  first  one  identifies  the  "center"  that  is  used  in  mapping  that  vertex 
onto  M  Thus  after  the  initial  vertex  0 .  the  nodes  7,...,1 2  become  centers  which  serve  to 
map  the  nodes  1 3,. ..,42  onto  M.  Then  the  process  continues  with  nodes  1 7, 1 8, 1 9, 23, 
24, 28, 29, 33, 34, 38, 39, 42  as  centers.  This  is  no  longer  shown  in  the  figure,  but ,  in 
practice,  we  have  always  continued  through  this  further  stage.  It  resute  in  the  earlier 
mentioned  total  of  1 1 4  triangles  involving  1 9  centers  and  hence  as  many  Jacobian 
evaluations 

Once  a  standarized  patch  ro(0  is  used  in  step  (5c),  a  suitable  divergence  check 
has  to  be  built  into  the  Gauss- Newton  process.  If  in  step  (5d)  this  check  is  triggered, 
then  the  corresponding  vertex  ^  of  2  is  flagged  as  unusable.  Such  unusable  vertices 
are  excluded  from  the  further  computation.  A  similar  procedure  may  be  followed  when 
in  step  (5b)  the  moying  frame  algorthm  fails  However,  in  the  latter  case  it  is  often 
advantageous  to  re-mtialize  the  moving  frame  algorithm  at  one  of  the  successfully 
computed  points  x  on  M  Of  course,  then  the  computed  basis  U(x)  has  to  be  used  as 
the  reference  matrix  T0. 

The  above  provisions  may  result  in  triangulations  that  cover  a  somewhat  irregular 
domain  on  M  Fortunately,  in  practice,  this  does  not  ocar  as  frequently  as  might  be 
expected,  provided  the  step  length  6  is  not  chosen  too  large.  A  Ornately,  t  is  natural  to 
consider  a  procedure  which  adapts  S  in  cases  of  failures  and  hence  which  produces 
irregular  tria regulations  on  M  Such  a  program  is  now  being  implemented. 


5  Numerical  Experiments 

In  mts  Section  we  present  resuls  of  some  numerical  experiments  with  the 
triangulabon  algorithm  The  method  produces  a  weath  of  numerical  data  which 
cannot  be  reproduced  in  the  limited  space  of  this  paper  At  the  same  tme,  tea 
challenging  problem  to  invent  instructive  graphical  representations  of  manifolds  of 
dimension  larger  than  2  As  a  consequence,  only  some  graphical  results  for 
two-dimensional  ma  nifolds  are  shown  here  t  is  hoped  that  other  experiments  with 
higher  dimensional  manifolds  can  be  presented  elsewhere 

An  Exothemr.  Reaction:  As  a  first  example  we  consider  a  simple  transport  model 
for  an  exothermic,  first-order  reaction-scheme  discussed  in  [4]  which,  in 
dmensionless  form,  leads  to  a  two-point  boundary  value  problem 

(Duy  ♦  ko0i-u)exp(-A(Uu)-',)  =  0,  u(0)=u(L)  =  0  (5.1) 

The  dmensionsless  parameters  and  A  involve  the  constant  concentration  and 
temperature  on  the  outside  of  the  system,  and  for  the  calculation  we  follow  [4]  and  set 
D*  1,  L  =  1,and  1^=  10^  The  standard  finite  difference  approximation  of  (51)  on  a 
uniform  mesh  Xj«  ih .  i  =  0,1,  ,n+1.  h  =  (n+1)"1  then  produces  a  nonlinear  equation 
of  the  form 

-«,.1  +  2xrxl+1  =  h-k0(|i-x,)exp(-A(1 +X,)"1) ,  i=1 ,n.  xo  =  xn+1  =  0  (5  2) 

For  p  =  1  there  is  a  simple  turning  in  A  near  A  =  22  which  was  calculated  with  the 
continuation  code  PI  TOON  and  n=1 0  This  point  was  then  used  to  initialize  the 


trangulation  algorithm  Here  --  as  in  the  subsequent  example  --  the  reference 
trangulation  in  was  generated  by  repeated  pivoting  by  reflection  from  the 
equilateral  triangle  (4  4).  The  stepsize  in  the  a/fine  mapping  (4.5)  was  6  =  0  4  Figure  2 
shows  the  resuls  of  this  tnangulabon  More  specifically,  the  intersection  of  the 
computed  smplicial  approximation  wth  the  (A,  ji,  x*-)  -space  is  shown,  where  xc  is 
the  computed  x-yalue  at  the  center  of  the  interval  The  coordinate  axes  are  slightly 
rotated  and  one  see  clearly  the  stars  taped  pattern  of  the  reference  tria  regulation  of 
Figure  1  as  it  was  mapped  onto  the  manifold  The  floor  of  the  valley  actually 
represents  a  foldline  as  can  be  seen  in  Figure  3,  where  the  same  surface  in  (A,  ji,  xj  - 
space  is  projected  onto  the  (A,  p) -parameter  plane,  t  shows  also  that  the  location  of 
the  turning  point  in  A  depends  approximately  linearly  on  ji 

A  Shallow  Crcular  Arch  As  the  second  example  we  consider  a  finite-element  model 
for  the  deformation  of  a  thin,  shallow,  crcular  arch  which  has  been  used  as  a  test  case 
by  many  authors  Jt  appears  to  go  back  to  A  C  Waker  [1 9)  and  we  employ  here  the 
same  formulation  as  in  [1 3]  In  particular,  in  a  (r,©)-polar  coordinate  system  with  the 
vertical  drection  as  the  r-axis.  the  unloaded  configiration  of  the  arch  is  represented  by 
the  crcular  segment  {(r,©).  r=  10, -©0s  ©s  60,  ©0  =  15°)  and,  for  pinned  ends, 
the  dimensionless  total  potential  energy  and  associated  boundary  conditions  are  given 
by 


+  ctj  (u")  -  ot^p 


\ 


d© 


u(  ©)=w(  ©)=  u"(  ©)=  0,  ©=  ±  © 


where  primes  denote  derivatives  with  respect  to  ©  For  the  asymmetrical  bad 


2  A 


p(9)  =  (1  -n)A ,  t  -e0  i  ©  <0,  and  p (©)  =  (1  +ul)A,  1 0  s  e  s  e0, 

involving  the  two  parameters  p.  and  A,  the  load  path  for  ^  =  0  has  a  bifurcation  point 
near  A  *  1.9  This  point  was  computed  wfth  PITCON  and  used  to  initialize  the 
triangulation  algorithm  The  stepsize  of  the  mapped  triangles  was  S  =  0.5.  The  results 
are  shown  in  Figure  4  More  specifically,  let  ^  denote  the  (dimensionless)  radial 
deformation  at  the  center,  then  the  figure  shows  the  intersection  of  the  manifold  with 
the  (xc,A,p)-space  projected  orthogonally  onto  the  (A,ji)-plane.  The  cusp  -  bifurcation  is 
clearly  visible  and  the  saddle  shape  of  the  surface  can  be  seen  even  belter  in  the 
slighty  rotated  Figure  5. 

The  problem  was  also  run  wih  the  load  function 

p(©)  =  A(l-4(n-©X©0r1i  if  maxf-e^-O^S©^  *  ©  *n, 

P(©)  =  A[1+4(pi-©X©0r1t  if  \x  *  ©  *  min(©o,n+0.25©o), 

considered  already  m  [3]  In  other  words,  the  load  is  a  piecewise- linear  hat  function 
which  has  the  value  A  at  ©  =  pi  and  is  zero  outside  the  interval  centered  at  p.  of  width 

0  5©q 

Figures  6  arid  7  give  results  with  this  load  obta  ined  at  two  initialization  points  More 
specifics  lly,  the  computations  for  these  figures  were  centered  at  the  limit  points  with 
respect  to  A  when  p  is  fixed  at  p=0or  p  =  0.05,  respectively  Once  again  these  limit 
points  were  computed  with  PITCON  The  fold  line  in  the  (p.A)-plane  has  the  shape 
shown  in  Figure  8  and  Figures  6ano  7  c  (early  s  how  seoments  of  this  fold  line  In 


particular,  Figure  7  contains  one  of  the  points  where  the  most  dangerous  load  occurs 
namely  at  about  \i  =  0.1 6©0. 


W 
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6  Outlook 


The  numerical  examples  indicate  that  the  tna  regulation  algorithm  works  very 
efficiently  even  around  singularities  Thus  -  as  intended  does  indeed  offer  a  new 
tool  for  denying  information  about  the  shape  and  features  of  the  mantfold  Of  course,  as 
mentioned  before,  we  are  able  to  present  here  only  some  graphical  information  and 
none  of  the  extensive  numerical  output  of  the  algorithm  This  output  is  available  as 
input  to  various  post-processes  for  extracting  further  information  Several  such 
processes  are  now  under  development  and  'will  be  desorbed  elsewhere  in  more 
detail  Here  we  present  only  a  brief  overview  of  some  of  the  possblities 

As  noted  earlier,  linear  interpolation  between  the  vertices  of  the  computed 
triangulation  defines  a  smnplicial  approximation  Ms  of  the  corresponding  part  of  the 
manfold  M  The  points  of  M2  can  be  used  to  compute  further  points  of  M  For  example, 
we  may  project  any  such  point  onto  M  by  applying  the  corrector  leration  (2.1 1). 
Alternately,  we  can  augment  the  system  (2.2)  with  p  appropriately  chosen  equations 
and  then  apply,  say,  a  chord  Newton  method  to  calculate  a  corresponding  point  on  M. 
This  approach  is  useful  when  points  on  M  wth  specfic  properties  are  desired.  For 
example,  we  may  be  interested  in  certain  target  points  where  the  parameters  have 
prescribed  values  In  that  case,  these  target  conditions  become  the  augmenting 
equations  and  we  may  start  the  ierative  process  from  a  point  on  where  the 
parameters  have  the  speefied  values  Augmenting  equations  are  also  essential  when 
we  are  interested  in  determining  the  specfic  location  of  certain  types  of  singularities 
For  the  computation  of  limit  points  a  comparison  of  various  such  augmentations  was 
given  in  [1 3],  and  for  higher  order  singularities  the  literature  on  suitable  augmentations 
is  fa  rly  extensive  In  our  present  setting,  the  "minimal”  augmentations  discussed  in 


[Ill,  and  thereafter  in  [10]),  appear  to  be  of  particular  interest 

For  any  given  functional  the  computed  data  allow  us  to  generate  contour  plots  on 
the  intersection  of  the  srnplicia!  approximation  with  various  subspaces.  For 
instance,  in  some  structural  problem  we  might  be  interested  in  seeing  lines,  where 
some  stress  component  is  constant,  plotted  in  dependence  of  certain  other  variables 
A  special  example  of  such  a  contoir  plot  involves  the  graphical  representation  of  lines 
of  foldpomts  As  the  Figures  of  the  previous  section  already  show,  our  triangulations 
provide  information  for  approximating  segments  of  such  foldlines  One  approach  for 
detecting  foldpomts  is  to  monitor  the  orientation  of  the  projection  of  the  tangent  basis 
onto  the  parameter  space  t  there  is  a  change  in  this  orientation  then  we  have  passed 
through  a  foldpoint,  but  the  converse  is  not  necessarily  true,  that  is,  not  every  foldpoint 
can  be  detected  this  way  The  orientation  is  characterized  by  the  dete.minant  of  the 
projected  basis  in  the  parameter  space  Thus,  if  we  plot  lines  of  constant  determinant 
values,  then  lines  of  zero  determinant  are  approximations  of  the  desired  foldlines.  Of 
coirse,  for  this  we  need  the  tangent  basts  at  each  vertex  of  the  Iriangulation  and  that 
increases  the  cost  of  the  overall  algorithm  However,  there  appear  to  be  other  possible 
techniques  for  approximating  foldlines  from  data  obtained  by  oir  triangulation 
algorithm  This  will  not  be  pusued  here 

Even  though  these  contour  plots  only  provide  lines  of  constant  values  on  the 
smplicial  approximation  rather  than  on  M  itself,  they  tend  to  offer  already  good 
insight  into  the  shape  of  the  ma  nifold.  Of  course,  as  discussed  earlier,  we  can  always 
call  on  various  local  corrector  methods  to  project  these  lines  onto  M  itself 

So  far  we  mentioned  only  the  need  for  appropriate  post-processing  techniques 


for  analyzing  the  output  of  our  trianguiation  method.  There  is  also  considerable  room 
for  mproYing  the  algorithm  itself  In  particular,  for  large  sparse  problems  the 
QR-factorization  may  be  computationally  expensive.  As  noted  earlier,  there  exist 
resuts  for  computing  sparse  bases  of  the  null  space  of  a  matrix  (cf  [5],  [6]).  The  rotation 
requr*  ’  for  the  moving  frame  algorithm  is  Hkely  to  destroy  this  sparsiy,  and  hence  the 
rotated  basis  should  not  be  stored  but  computed  only  as  needed.  In  the  case  of  low 
drnensional  manifolds,  this  is  not  very  costly  as  long  as  the  computation  of  the  original 
basis  of  the  tangent  space  takes  account  of  the  sparsity  structure  of  the  Jacobian 
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